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^r^ Abstract. Particles that are immersed in a fluid exchange momentum via the 

c/3 ■ fluid, hence their Brownian motion is correlated. By means of multiparticle-collision 



dynamics simulations we study the interactions between two colloidal beads in a 
sheared fluid suspension. Recently, this topic has been addressed in experiments 
on colloidal particles trapped by optical tweezers in a microfluidic device [PRL 103, 
^fj I 230602 (2009)] and theoretically by means of a Langcvin model [Eur. Phys. J E 33, 

Ch ■ 3f 3 (2010)]. Although we neglect the rotational degrees of freedom of the colloids, and 

employ a very simple coupling between the colloids and the flow field, we can reproduce 
the experimental data and partly explain why it differs from theory 
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1. Introduction 

A mesoscopic particle that is suspended in a liquid moves stochastically due to collisions 
with the particles of the liquid - it performs Brownian motion. When two particles are 
suspended in a liquid, they transfer momentum from and to each other via the liquid, 
hence their Brownian motion becomes correlated. This coupling of stochastic motion 
via hydrodynamic interactions is important in biological systems as well as microfiuidic 
devices. 

Recently Zimmermann and co-workers presented a series of articles on the 
correlations that are induced between two trapped colloidal particles by shear flow. They 
addressed the problem experimentally [1] and theoretically [2J [3]. The experimental 
system consisted of uncharged polystyrene beads suspended in water in a microfiuidic 
device that allowed to create a linear shear profile. The particles were trapped by optical 
tweezers and their positions were recorded by with high speed camera. We will refer to 
the results of the experiments in section 14.21 In the theoretical work a Langevin model 
was solved. We briefly introduce the model and its solutions in section [21 

Here we present a computer simulation study on the same problem. In section 
|3] we describe the simulation method, and in section H] we compare our results to the 
experimental and theoretical results. 

2. Equations of motion 

We briefly review the work of Zimmermann and co-workers [2, [3]: The systems 
condsidered here are that of one and of two hard spherical particles trapped in harmonic 
potential wells and suspended in a shear flow with a linear velocity profile. Here and in 
the following the Brownian spheres have the same mass Mi = M and the same effective 
radii Qi = g with i = 1,2. Their position vectors are denoted by r^ = [x^y^zi). The 
isotropic potential wells Ui (rj) = | (t*j — qi) are located at q x = (0, 0, 0) (single bead) 
and qi = (±D/2, 0, 0) (pair of beads). Here k is the strength and D is the separation of 
the potentials. The potentials give rise to forces given by Ff = — VC/j. As the colloids 
are surrounded by a fluid, they experience a force F f due to the friction with the fluid 
molecules. For small Reynolds numbers Re this force is proportional to the relative 
velocity of the colloid and the local fluid field, F? = —( (u — ri), where ( = 6nr]g is the 
Stokes friction coefficient and r\ is the fluid viscosity. The linear velocity profile of the 
sheared suspension is given by uir) = ^ze. x with the shear rate 7. We are interested 
in the case of overdamped motion, i.e. the inertial part Mri is negligible, and one can 
write down the Langevin equation 

r { = u{n) + HijFj + Fi (1) 
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Here, r 12 := t*i — r 2 is the distance vector between the two particles and its absolut 
value r 12 = \r\i\. The last contribution in eq. ([Q) is the stochastic force acting on 
particle %. For a fluid, in which orthogonal components of velocity fluctuations of the 
fluid molecules are uncorrelated, the following characteristics of F* can be assumed: 

{Fm) = °' (3) 

(F l s (t)F^t')) = 2k B TH lJ 5(t-f). 

One can rewrite eq. JT} in a more compact form (cf. [2]) by defining the vectors 
R := (7*1, r 2 ) and Q := (gi, g 2 ) as well as the 6x6 matrices 

and U(R) := UR, where U13 = U^ — 7 and all other Um = 0, so that it now reads as 

R=UR + kH(Q-R) + F 1 (5) 

in which F satisfies © with Ff -> F = (F?,F%) and H tj -> H. The coupled 
equations of motion (jSJ) have been solved by Zimmermann and co-workers (2j [3]. Their 
solutions were given in terms of the particle fluctuations around the potential minima 
R(t) = R(t) - Q as 

R(t) = e~ tM R(0) + f dt'e^'^F, (6) 

Jo 

where M := kH — U . They form the basis of the analytic calculation of the correlation 

functions (R(0)R(t)), where the brackets (. . .) denote an ensemble average. 

3. Simulation details 

In order to simulate colloidal particles embedded in a sheared fluid environment, we 
used a combination of a simple molecular dynamics (MD) simulation and multiparticle- 
collision dynamics (MPCD), a mesoscopic solvent model to account for hydrodynamic 
interactions [3]. 

3.1. Multiparticle- collision dynamics 

The MPCD method we employ has been developed to solve the equations of 
hydrodynamics in a fluctuating solvent. N point particles of the same mass m are used 
for the transport of momentum through the system, while satisfying the conservation 
laws of mass, energy and momentum locally. The algorithm consists of two steps, namely 
free streaming interrupted by multiparticle collisions. In the streaming step, all fluid 
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particles are propagated ballistically with their velocities v k , i.e. the particle positions 
are updated according to a time increment h as 

x k (t + h) = x k (t) + hv k (t). (7) 

After the time step h, the N fluid particles are sorted into a lattice of cubic cells of size 
ax ax a, so that on average n particles are within each collision cell. Then, the particle 
velocities are rotated around the center of mass velocity v in this cell, 

v k (t + h) = v(t) + fi(a) [v k (t) - v(t)} , (8) 

where Q(a) is a rotation matrix corresponding to a fixed angle a, which is generated 
randomly for each cell. Before a collision step is carried out, the collision cell grid 
is shifted by a randomly chosen vector with components taken from the interval 
[—a/2, a/2], to ensure Galilean invariance [5]. As the transport coefficients like the 
dynamic viscosity n can be expressed as functions of the simulation parameters in an 
analytical form j6j[7J[8], one can tune their values to fullfill the conditions of overdamping 
and low Re given in sec. [2j 

The coupling between the colloidal particles and the fluid particles is accomplished in 
the simplest manner, which means that a point particle with instantaneous velocity 
V and mass M takes part in the collision step eq. flH]) within its cell with v = 
(MV + m E" v)/(M + nm). 

In order to shear the fluid in our system, we confined the simulation box of size LxLxL 
using two flat walls in z direction at ±L/2. Those were moved in x direction with 
velocities t4 = iijL/2. To reduce fluid slip at the walls we used a bounce-back rule, 
i.e. if a fluid particle hits a wall during the streaming step, its velocity is inverted in the 
rest frame of the wall (v 1 — > —v'), where «':=«- «; . In combination, we used an 
algorithm prosposed in [9]. Here, N pp pseudo-particles are inserted into the cells, which 
are cut by the flat walls as a consequence of the random shift, so that in the partially 
filled cell the average number density is restored, n + N pp = n. 

3.2. Molecular dynamics 

In the MD part of our simulations, the positions of the colloidal particles were computed 
according to the Velocity Verlet algorithm with time step At several times between two 
consecutive MPCD steps. The number of position updates is given as h/At. In the case 
of a pair of beads, we did not account for collisions between those. Hence we would 
expect deviations from the behaviour derived in sec. [2] for small distances between the 
potential wells. As the "colloidal" beads also took part in the collision step as point 
particles, there is no physical volume associated with them, and therefore no geometrical 
radius can be defined. 

For our simulations we chose the following parameters: the temperature was set to 
ksT = 1 and was kept constant on average using the thermostat proposed in ref. [TO] . 
The collision cell size a and the solvent particle mass m were set to unity as well. The 
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simulation box had dimensions 30 a x 30 a x 30 a. The mass of the colloids was set 
to M = 2m. The rotation angle had the value a = 130°. The two time steps of 
our hybrid dynamics simulations were At = O.OOOlr and h = O.Olr, both in units of 
r := ^Jma? /ksT. The collision cells were occupied on average by n = 10 point particles. 
With those parameters, the dynamic viscosity was calculated as n « 82.2 ^/kBTm/a 2 . 
The value for the strength of the potential wells was k = 10 ksT/a 2 and the rate with 
which the fluid is sheared was set to j — 0.04 r _1 . With those values, a Reynolds num- 
ber of O(Re) ~ 0.1 <C 2300 results. Furthermore, the ratio (/m, which is a measure for 
the damping in the system, is of the order of O (C/m) ~ 10 2 3> 1. Hence, the trapped 
particles perform an overdamped motion embedded in a laminar fluid flow. 

4. Results 

In the following we will present our results for both the auto-correlations (AC) and the 
cross-correlations (CC) in the random displacements of Brownian particles embedded 
in a liquid. The simulations were done in a quiescent as well as in a sheared fluid. We 
focus on fluctuations parallel (x) and perpendicular (z) to the direction of the shear 
flow. The correlation data was measured as the average of 10 5 — 10 6 time series. Error 
bars are given by the standard deviation of the data sets. 

4-1. Single sphere 

Given the defining properties of the stochastic forces (eq. (jBJ)), the Langevin equation 
eq. ([1]) for a single particle embedded in a fluid without flow leads to an isotropic AC 
function, which decays exponentially with a relaxation time r p = C,/k 

(x(o)x(t)) = (mm) = ^f e ~ t/Tp - ( 9 ) 

Using this equation, we determined the relaxation time from our simulation data as 
r p = 23.7 r. As the fluctuations of orthogonal fluid velocity components do not couple 
for 7 = 0, the CCs vanish, (x(0)z(t)) = (x(t)z(0)) = 0. 

The situation changes if the colloidal bead is exposed to a sheared fluid. Then the AC 
in the shear direction is modified to 

k B T 



(x(O)x(t)) 



k 



Wi 2 / t 



e~ t/Tp , (10) 



2 

while (z(0)z(t)) equals the one of eq. ([9]). Here, eq. ( TTUj) is given in terms of the Weis- 
senberg number Wi := 7T P = 0.95. Both correlation functions are shown in figure [TJ 

Now, if there is a finite shear rate, the CCs are not zero anymore and fluctuations 
of orthogonal components couple under the influence of the linear velocity profile of the 
surrounding fluid. This can be visualized by the distribution of the particle positions 
in the shear plane, as it is depicted in figure [2j Here the distribution, which would be 
spherical in the case of a quiescent fluid, assumes an elliptical shape. As it was shown 



HS under shear flow 



o < x(0)x(t) > 
n <z(0)z(t)> 
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Figure 1. Auto-correlation of particle position for a single particle in the direction 
of shear flow (simulation: black circles, theory: green line) and perpendicular to it 
(simulation: red squares, theory: blue line). The arrow denotes the size of the error 
bars. 



in ref. [3], the characteristic dimensions of the ellipse, as the ratio between the principle 
axes w/l and the inclination angle 0, are connected to Wi via 

1 



tan i 



w/l 



4 + Wi z - Wi 




1/2 



;n) 



Wi 2 + Wi, 

From the data in figure [2} we found = 32.3 ° and w/l = 0.622. Using equations ( fill , 
those lead to the Weissenberg numbers Wi = 0.95 and Wi = 0.98, respectively, both in 
good agreement with the value calculated above. Also, in the simulation the CCs do 
not vanish any more. They show a behaviour which satisfies the equations 



(x(t)z(0)) 
(x(0)z(t)} 
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Both functions are plotted in figureS There, the initial increase of (x(t)z(0)) reflects the 
fact, that finite fluctuations z(0) 7^ are carried away by the shear flow in x-direction 
before the initial displacement starts to relax. The analytical predictions are in excellent 
agreement with the simulation data (no fit parameters). 



4-2. Two hard spheres 

In the case of two trapped particles, not only correlations between fluctuations in 
different directions of each colloid are of interest, but also inter-particle correlations 
like (xi(t)x2{0)) . Those are functions of the distance between the potential minima 
D = \Q12\1 where q\2 '■= q\ — <?2 is the connection vector of the two potential wells. In 
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Figure 2. Particle position distribution in the shear plane for a single particle in 
sheared flow. The box denotes the size of a collision cell ax a. 



o < z(0)x(t) > 
d < x(0)z(t) > 




time [Tp] 

Figure 3. Cross-correlations of particle position for a single particle in shear flow. 
The arrow denotes the size of the error bars. 



order to compare our results with analytical expressions, one has to be in a regime where 

the hydro dynamic radii g of the particles are small compared to D. To obtain a value 

for g, we used the definitions of the relaxation time r p = (/k and the friction coefficient, 

from which it follows that g = Tpk/Gnr]. In our simulations the hydrodynamic radius 

roughly equals 0.15 a. 

In this section we will follow the nomenclature of ref . [2] and introduce the four relaxation 

rates: 

1-2// 



Ai 



A 2 
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with the parameter /i := 3g/AD. To illustrate the meaning of the Aj, let us consider a 
situation in which the two particles, after experiencing independent stochastic kicks, are 
pulled back into the minima of their potential wells, respectively. This relaxation process 
can be decomposed into two modes for each spacial direction, namely a parallel and an 
anti-parallel translation of the particles. As the hydrodynamic interaction between the 
spheres depends on their relative motion, restoring is accelerated in the parallel and 
damped in the anti-parallel case. Here, the modes for the two directions perpendicular 
to qi2 coincide, so that only the above relaxation rates are left over. 
A previous experimental study [IT] with two trapped Brownian particles in a quiescent 
fluid showed anti-correlations between random fluctuations along same directions. These 
anti-correlations are a function of \i. In addition, the functions are modified under the 
influence of a linear shear profile in a similar manner as the AC function for fluctuations 
in shear direction in the one-particle case ( JTOj) . Including this shear flow correction, 
which is proportional to Wi 2 , {xi(0)x 2 (t)) is given by 

(xi(0)x 2 (t)) =\{e- X ^ - e" A3i ) 

Wi 2 / - (1 + /i) e~ Xlt e~ X2t 

7 V 6« 2 + 7« + 2 + 2 + 3/- 

Wi 2 / - (1 + fi) e~ X:it e~ Mt 



+ 4/i V 6/i 2 + 7/i + 2 + 2 + 3/J 



4« V 6« 2 -7« + 2 2-3/x, 

The simulation data for the CC functions and eq. (fT4l) are plotted in figure |4] for an 
unsheared system (Wi = 0) with two different D and a sheared (non-zero Wi) system. 
Again, the agreement with the analytical predictions is very good. We note that not only 
the depths of the minima differ but also the position of the strongest anti-correlation 
changes when the particles are embedded in a sheared environment. 
In the data for Wi = and D = 0.5a (black circles in figure H]) there is an additional 
decrease of the CC function at long times, which is not covered by the error bars. 
A similar behaviour can be found in the data of ref. [11] for small distances between 
the laser potentials. The origin and the physical relevance of this effect need further 
investigation. 

Finally, we will briefly discuss inter-particle correlations between perpendicular 
fluctuation directions, namely (xi(0)z2(t)} and (zi(0)x2{t)) ■ The behaviour of the latter 
can be understood as a combination of the inter-particle anti-correlation of fluctuations 
in shear direction and the fact, that a fluctuation in ^-direction of a single particle is 
followed by a positive motion in x-direction (cf. equation ( TT2l) ) 

~ antiCC ~ ~ / ~ \ shear flow _ 

x 2 — V xi = xi Ui) — > zu 
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Figure 4. Cross-correlations of fluctuations parallel to the shear direction, 
(ii(0)x2(t)), between two beads for unshcarcd (simulation: black circles, blue 
diamonds, theory: green and orange line) and sheared system (simulation: red squares, 
theory: blue line). 
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o <z 1 (0)x 2 (t)>, D = 1.5a 



time [Tp] 

Figure 5. Cross-correlations of fluctuations in perpendicular directions between two 
beads for D = 0.5a (black circles, green line and red squares, indigo line) and D = 1.5a 
(blue diamonds and orange line). 



which produces an anti-correlation in (zi(0)x2{t)). The analytical expressions for the 
two inter-particle CC functions are given by 



. . , . ,. Wi / e" A2t 
( Xl (0)z 2 (t))= — 



4 V 2 + 3/i 2-3fi 

Wi , w 
&(0)5*(0>=— (-■ 
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These expressions are compared to the simulation data in figure [5j Again, they agree 
within the errorbars. 
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4-3. Comparison with experiment 

In this section we compare our results to the experimental data from ref. [3]. The 
green stars in fig. [6] are the experimental values for the single particle cross correlation 
functions (x(0)z(t)) (left panel) and (z(0)x(t)) (right panel). There are clear differences 
with respect to the simulation results in fig. El (x(0)z(t)) shows a pronounced minimum, 
while the simulated curve drops monotonically (red squares in fig. [3]), and the decay of 
the experimental (z(0)x(t)) is much faster than of the simulation data (black circles in 

fig. ED. 

In contrast to our simulation (in which there was only one particle), in the experi- 
ment the single-particle CCs were obtained using a system that contained two particles. 
Hence, the single-particle correlations of perpendicular fluctuations are influenced by 
presence of the second particle. For small potential well distances this can lead to de- 
viations from the analytical description given in j2], which is only valid for small /x, i.e. 
large distances D. In the PhD thesis of A. Ziehl [12] it is suggested that these deviations 
are possibly the origin of the additional minimum in (x(0)z(t)). To test this hypothesis, 
we also analysed the CCs for reasonably large \i with a second particle present. The 
simulation data is plotted in figure EJ Here, the functions were renormalized to account 
for the different Weissenberg number of Wi exp ~ 0.62. Our data shows the development 
of a minimum with increasing ji which supports the hypothesis. However, the minimum 
is much less pronounced than the one in the experimental data. Hence there might be 
other sources for this effect. 

Fig. [7] shows the experimental data for the two particle correlation functions. There 
is a clear difference in time-scales between simulation and experiment (left panel). This 
is suspected to be due to deviations in the prodcution process of the microfluidic device, 
which cause changes in the trapping potential [13] . To compare our results of the 
inter-particle CC functions, we renormalized the time scale to match the experimental 
relaxation time r oxp , which is given by the position of the minimum of (zi(0)x2{t)) and 
which differs from the relaxation time used before in ref. [3] . We again took into account 
the different Wi. After these transformations have been applied, the data sets show very 
good agreement (right panel). 

5. Summary 

We have presented a computer simulation study of the hydrodynamic interactions 
between two colloids in a sheared fluid. In particular, we have tested the hypothesis 
that the recently observed non-monotonic behaviour in the cross-correlation function of 
the position of a single particle is due to the presence of a second particle. We found 
that a second particle produces this effect, but it might not be the full explanation. 
Another source could be the rotational motion of the colloids, which we neglected in 
the simulation. 
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Figure 6. Singlc-particlc cross correlation functions (x(0)z(t)) (left panel) and 
(i(i)z(O)) (right panel). Comparison of simulation data for /i = 0.0038 (black circles) 
and jjl = 0.228 (red squares) with experimental results (green crosses). 
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Figure 7. Inter-particle cross correlation functions. Comparison of simulation data 
(black circles and red squares) with experimental results (green crosses and violet 
stars) . Left panel shows data with time rescaled by r p . Right panel shows results with 
time rescaled by r Gxp , given by minimum of experimental (zi(0)x2(t)) , and correlation 
strength rescaled with respect to different Weissenberg numbers. 
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